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Forward time step integrators are splitting algorithms with only positive splitting coefficients. 
When used in solving physical evolution equations, these positive coefficients correspond to positive 
time steps. Forward algorithms are essential for solving time-irreversible equations that cannot be 
evolved using backward time steps. However, forward integrators are also better in solving time- 
reversible equations of classical dynamics by tracking as closely as possible the physical trajectory. 
This work compares in detail various forward and non-forward fourth-order integrators using three, 
fourth, five and six force evaluations. In the case of solving the 2D Kepler orbit, all non-forward 
integrators are optimized by simply minimizing the size of their backward time steps 



I. INTRODUCTION 



Many physical evolution equations are of the form 



-g^ = {T+V)w, (1.1) 



where T and V are noncommuting operators. Important examples include the imaginary time Schrodinger equation 

dr 4 



9^ (V-y)V^ (1.2) 



and the Fokker-Planck equation 

^p(x, t) = 1 V2p(x, <) - V • [v(x)p(x, t)]. (1.3) 

Because the diffusion kernel cannot be evolved backward in time, both of these are time-irreversible evolution 
equations. Aside from these obvious examples of p.ip . any pair of equations of the form 

can also be casted into the form (|l.ip . This is because the evolution of a general function W{p,q) (including q and 
p themselves) can be formulated as 

dW _ dW dq dW dp 
dt 9q dt dp dt 

[^{p).^+F{q)-^)w, (1.5) 



from which one can identify 



9q dp 



T^ = v(p)-^ and F = F(q).— . (1.6) 



Classical Hamiltonian dynamics corresponds to 



v(p) = P (1.7) 



and the resulting evolution (|1.4p is time-reversible. 

The generic evolution equation (|l.ip can be solved iteratively 



w{t + e) =e''''^+^'^w{t) (1.8) 
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by approximating the short time evolution operator e^'"^+^^ to any order in e via 

N 

e-(T+v) ^■Qgt.sTg..sy^ (1.9) 
1=1 

assuming that the effect of e^"^ and e^^ can be computed exactly. The set of coefhcients {ti,Vi\ are determined 
by the order condition. For classical dynamics (|1.6p . every factorizations of the form (|1.9p produces a symplectic 
integratori^i^iiii^i^i^i^ii^ which is an ordered sequence of alternating displacements of p and q preserving Poincare 
invariants. For periodic motion, their energy errors are bounded and periodic, in contrast to explicit Runge-Kutta type 
algorithms whose energy error grows linearly with the number of period o^^i^^ . However, for solving time-irreversible 
evolution equations such as (|1.2[1 or (|1.3p with T = V^/2, the Green's function 

G(r',r;t,£)cxe-('^'-'")'/(2t.e) (l^^O) 

is the diffusion kernel only if ti is positive. If ti were negative, then the kernel is unbound and there is no way of 
simulating the diffusion process backward in time. 

Historically, symplectic integrators were developed extensively for use in classical dynamicai^i^ii^i^i^i^iii^ii^iiS. Since 
classical dynamics is time-reversible, there was no impetus for demanding that all ti be positive. Moreover, Shengi^ 
and Suzuki^^ have proved that all factorizations of the form ()1.9|) beyond second order must contain some negative 
coefficients in the set {ti^Vi}. Goldman and Kaper^^ later proved that for factorizations of the form (|1.9p beyond 
second order, both operators must have at least one negative coefficient. Thus all conventional splitting schemes beyond 
second order must contain some negative coefficients and none can be used to solve time- irreversible problems. Because 
of the Sheng-Suzuki theorem, it is also difficult to see how one can devise all positive coefficients, forward time-step 
algorithms. 

The operator product p.9p has the general Campbell-Baker-Hausdorff expansion, 

n e*-^^e"'^^ = expe^CTT + eyV + eervlT, V] + e^eTTv[T, [T, V]] + e^evTv[V, [T, ¥]] + ■■ ^ (1.11) 
where all the error coefficients er, erv, ^vtv, etc., are calculable functions of {ti,Vi}, in particular, 

N N 

bt — ''^^ti and ey^^^w^. (1-12) 

i=l i=l 

In order for the product to be consistent with the original evolution operator, the coefficients {ti,Vi} must satisfy 
the above constraints with ct = 1 and ey = 1. Forcing the remaining error coefficients to vanish results in order 
conditions that {ti,Vi} must satisfy. It is easy to force ctv = 0- Any left-right symmetric product will do. For 
example, 

T2{e) = e^^ V^e^^^ ^ ^e{T+V -:i^e^[T.\T,V]] + ^e'[V.\V.T]] + -) -^3^ 

produces the following second order symplectic algorithm according to (|1.6p : 

, 1 Po 
qi = qo + — 
2 m 

Pi = po+eF(qi) (1.14) 

^1 Pi 
q2 = qi + — 
2 m 

where the last numbered variables are the updated variables. Thus any symmetric splitting with gt = ey = 1 will 
result in at least a second order algorithm. The surprise is that, as first shown by Shengi^, beyond second order a 
general sum of products of the form (|l.lip is incompatible with having positive coefficients. More specifically, Suzuki^^ 
shown that the two error coefficients cttv and evTV cannot both be forced to zero for positive coefficients {ti,Vi}. 
Since Takahashi and Imada^^ have shown that [V, [T,V]] = \W{r)\'^ is a local potential function when solving the 
imaginary time Schrodinger equation, Suzuki suggested^^ that this error commutator be kept and ways be found to 
efiminate [T, [T,V]]. 

Following up on Suzuki's suggestion, this author derived three simple fourth-order forward algorithms'^ in 1997 and 
demonstrated their efficiency in solving Kepler's orbit. Interestingly, it was found that classically [V, [T, V]] produces a 
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force (also with potential |Vl^(r)p) first derived by Ruth-'^ via canonical transformations. The two forward schemes A 
and B derived in Ref.^^ were also known to SuzukiJ^ based on McLachan's result?^ on slightly perturbed Hamiltonians. 
However, Suzuki did not implement them to do any calculation. 

Since 1997, fourth-order forward algorithms have been widely applied to time-irreversible systems such as the 
Fokker-Planck equation in deriving the first fourth-order Langevin algorithm^^, the Kramers equation^^ for describing 
stochastic dynamics, the Diffusion Monte Carlo algorith m^^'^'^ for solving quantum many-body ground states, the 
grid based imaginary time Schrodinger^^ equation in doing density functional calculations, Path-Integral-Monte Carlo 
method a^^i^^'^^i^^'^^i'^° for computing the quantum trace at finite temperature, short time evolved wave functions^ 
for doing variational quantum many-body calculations, the Gross-Pitaevskii equation for describing a Bose-Einstein 
condensate in a rotating anisotropic trap^ and electrons in a magnetic field confined by quantum dots^^ and rings24. 
These forward fourth-order algorithms are far more accurate than second order algorithms and allow very large step 
sizes to be used. 

While forward algorithms are indispensable for solving time- irreversible equations, from their inceptioniS, they 
have also been shown to be efficient in solving time-reversible equations. In comparison with explicit Runge-Kutta 
algorithms and conventional non- forward symplectic integrators, forward algorithms have been shown to be superior in 
solving the Kepler problem.^^'^^''^^, gravitational few-body problems^^, the real time Schrodinger equation^^, the real 
time Schrodinger equation in a laser field^^i^i^ and specially the radial Schrodinger equation^^. The incorporation 
of the commutator [V, [T, V]] in forward algorithms has also inspired a new class of higher order gradient symplectic 
integrators^2i43^^_ These new algorithms, through not forward beyond fourth-order, have less backward steps at 
higher orders. 

The reason why forward integrators are also better in solving time-reversible, classical dynamics problems is not 
well understood. From the perspective of the operator product approximation (|1.9p , as long as the second-order error 
terms are zero, any symmetric factorization scheme will be fourth-order; it should not matters whether these error 
terms are forced to zero with positive or negative coefficients. However, just as in the discussion of time-reversible 
and time-irreversible algorithms, one must move beyond the purely algebraic discussion of factorization schemes to 
examine how the resulting algorithms produce the solution of any particular equation. This work uncovers a crucial 
difference between forward and non-forward schemes when they are implemented as integrators for solving classical 
dynamics problems. For classical dynamics, if the trajectory is the exact solution to (|1.4p . then the force is evaluated 
only along the trajectory. As will be shown below, non-forward algorithms, when compared to forward algorithms at 
the same finite step size e, evaluate the force at intermediate positions far from the trajectory. As one reduces the 
size of the negative time steps, one also reduces the distance of these force evaluation points from the trajectory. In 
all cases examined, non-forward algorithms are improved by simply reducing the size of their backward time steps, 
allowing the force to be evaluated more closely along the trajectory. One can argue that these intermediate force 
evaluation points are not the trajectory outputs of the integrator and their placements are not required to be on 
the trajectory. That is correct. However, since the exact trajectory is determined only by forces evaluated on the 
trajectory, any unnecessary force evaluation off the trajectory is just "wasteful". One then must reduce the time 
step to bring the force evaluation points closer to the trajectory. This is exactly what is observed in the following 
comparisons. To achieve the same accuracy, non-forward integrators must use smaller time steps. The question here 
is not about correctness; it is about efficiency. In the following sections we will compare in detail both types of 
integrators with three, four, five and six force evaluations. 



II. INTEGRATORS WITH THREE-FORCE EVALUATIONS 



We will begin by comparing the efficiency of integrators in solving the 2D Kepler problem defined by the Hamiltonian 

H=Ip'-^. (2.1) 



(Henceforth, we will always normalize the Hamiltonian with kinetic energy p^/2 and m — 1.) This problem is an 
excellent benchmark because one can gauge an integrator's performance not just by its energy error but also by its 
orbital precession erro r^^i'^^'^^ . The latter is a more direct measure of the accuracy of the computed orbit. 

There are basically three fourth-order integrators that required only three- force evaluations: 1) the non- forward 
Forest-Ruth (FR) integratoriMi^, 

TpRie) = T2iaie)T2iaoe)T2{aie) (2.2) 



where 



1 2i/3 
"1 = 2 „ 21/3 "° " ~ 2-2i/3 ^ "^-^^ ■ (2-^) 
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2) The forward integrator AiSiii^: 



with V defined by 



corresponding to an effective force 



T^ie) = ei^^e^^^ei^^e^^^e^-^, (2.4) 



V = V+^e^[V,[T,V]] (2.5) 



F(q)=F(q) + -i£2V(|F(q)p). (2.6) 
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Transcribing each operator in (|2.4p yields the integrator 



Pi 


= Po H 


h^eF(qo) 


qi 


= qo H 


1 

-^epi 


P2 


= PH 


h^eF(qi) 


q2 


= qiH 


1 

-2£P2 


P3 


= P2 H 


h^eF(q2), 



(2.7) 



Starting with initial values po and qg , the updated variables are p = Pa and q = q2 ■ Algorithm A only requires two 
evaluations of the force and one evaluation of the force gradient. Recently, Omelyan^ has suggested that the effective 
force (|2.6p can be evaluated by extrapolation 

F(q)=F(q+le2F(q))+0(£4). (2.8) 

The resulting integrator remains angular momentum and phase- volume conserving and is nearly indistinguishable from 
a fully symplectic integrator when solving the Kepler problem. We will denote this use of an extrapolated effective 
force (j2.8p as algorithm A'. Algorithm A' only requires three force evaluations. 3) The Runge-Kutta-Nystrom (RKN) 
integrator—. If one expresses p = Pa and q = q2 directly in terms of po and qo, then (|2.7p reduces to 

q = qo + epo + ^£^ (Fo + 2Fi) 
6 

P = Po + ^£(Fo + 4Fi+F2). (2.9) 
6 

This is the form of the RKN algorithm, with Fq = F(qo), F2 = F(q2), but with Fi = F(qi). The conventional RKN 
algorithm is defined by Fq = F(qo), Fi = F(q'j) and F2 = F(q2), where 

q'l = qo + |po + ^(|)'Fo 

q^ = qo + epo + ^£'F(q'i) (2.10) 

are the estimated midpoint and final position respectively. 

Fig[T] shows the fourth-order energy error coefficients of these four algorithms at step size e — P/5000, where P is 
the period of an highly eccentric orbit with initial values q = (10, 0), p = (0, 1/10) and eccentricity e = 0.9. The error 
coefficient is obtained by dividing the energy error by at smaller and smaller e until a convergent curve emerges 
independent of s. The curve is further normalized by the initial energy. Thus each algorithm has a characteristic error 
coefficient, its error "fingerprint", in solving the Kepler orbit. For symplectic algorithms, this convergence is already 
set in when e « P/1000. For the RKN algorithm, the energy error curve after the mid period keep on lowering with 
decreasing step size, showing no sign of convergence at finite e. The error only spikes at mid period near the pericenter 
point. Non-symplectic integrator such RKN are characterized by an irreversible increase in the energy error after each 
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period. The error spike of the FR algorithm is nearly ten times as large as that of algorithm A. The extrapolated 
gradient algorithm A' closely matches that of A. 

Since all symplectic integrators have periodic energy errors, the energy error is not the most critical benchmark. 
The orbital precession error, as measured by the rotation of Laplace-Runge-Lenz (LRL) vector'^^j^di, is more dis- 
criminating. Figl2] shows the rotation angle of the LRL vector along the trajectory of the particle. The corresponding 
error coefficient is again extracted by dividing by e*. In this case, the RKN integrator shows the same convergence 
as the symplectic integrators. Thus the precession error coefficient is well-defined and irreversible for all algorithms. 
The FR integrator's error is three times as large as that of RKN and 10 times as large as integrator A and A'. When 
the force gradient is extrapolated, the energy error remains periodic, but the precession error can differs substantially. 
Here A"s error is smaller, but in other cases it may not be. We will revisit this point later. 

In order to understand the poor performance of the FR integrator, we track its approach toward the pericenter point 
3 at two time steps earlier at point 1, as shown in FigO The time step used here is £ = P/400. The FR integrator 
consists of three applications of T2{e), resulting in three overlapping triangles. The first application of 72(aie) begins 
at position 1, evaluates the force at Fl, and lands at la. This is the triangle 1-Fl-la. The application of the backward 
substep 72(aoe) begins at la, evaluates the force at F2 and brings trajectory back past the starting point to 16. This 
is the backward triangle la-F2-lb. The final application of 72(aie) begins at lb, evaluates the force at F3 and lands 
the trajectory at position 2. This is the final triangle lb-F3-2. Starting at position 2, the algorithm repeats its three 
overlapping triangles and zigzags its way to point 3. It is remarkable that FR can achieve fourth-order accuracy by 
such a tremendous zigzagging. Notice that as the FR algorithm tries to turn the "corner" near the pericenter 3, all of 
its force evaluation points are far off the trajectory. In FigUJ the positions where each algorithm calculates the force 
are plotted. The backward loops executed by FR far off the trajectory is conspicuous. This is the fundamental reason 
by all non-forward algorithms perform poorly. By comparison, forward algorithm A always evaluate the force and the 
force-gradient close to the exact trajectory. Since RKN is similar in form to A, it strays from the exact trajectory 
only near the pericenter point 3. 

III. INTEGRATORS WITH FOUR-FORCE EVALUATIONS 

Because the error of the FR integrator is uncomfortably large, there is an ongoing effort to construct better 
non-forward algorithms by use of more force evaluations. A non-forward fourth-order algorithm can be obtain by 
generalizing (|2.2p to 

N 



%{e) = Y[T2ia,e), (3.1) 



1=1 

provided that the coefficients are left-right symmetric satisfyin g*^'^'"^^ 

N N 

^ fli = 1 and al = 0. (3.2) 

i=l i=l 

Unfortunately, for = 4, there are no real solutions to the above equations. We will examine algorithms of the 
general form 

T/v/i = . . . exp(rto7') exp(£wiT^) exp(£iiT) exp(£W2^) exp(£i2r), (3.3) 

previously studied by McLachlan's^. Since the algorithm is left-right symmetric, only operators from the center to 
the right are indicated. For a fourth-order integrator, the order condition requires that 

i;i = i-W2, t2^^-Uivl, to = l-2(ii +t2), (3.4) 



1 / /9ti -A±2w 



w = ^i-l2h+%tl V2 = -^\IT\^ "^' ""^ j (3.5) 

and that the free parameter ti < 0. There are four solution branches for V2- The choice of 

121 , 

= ^,(12 -Vm)^ -0.299 
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1 / I9ti -4 + 2w , , , 

-2 = 7 i + (3-6) 



with 

4 V V 3ii 

reproduces McLachlan's*^ recommended algorithm. By simply reducing the size of the negative time steps ti, we 
obtain results as shown in FigOand FiglHl The choice of ti = —1/24 yielded simple analytical coefficients 

1 + V5 11- V5 

V2 = — - — and t2 = — — — . (3.7) 

Also shown are results of forward integrator G^: 

Tc^... cm\sV) exp(i£T) exp(^£F) exp(i£T), (3.8) 



where V is as defined in (|2.5[) with the same interpretation (|2.6p . Algorithm C uses three force and one force gradient 
evaluations. The force gradient can again be extrapolated by another force evaluation. The results are similar and 
will be omitted in this comparison. To understand the poor performance of non-forward algorithms, we again plot 
their force evaluation points in Figl?! Starting at position 2, McLachlan's algorithm evaluates the force at Fl, back 
tracks and evaluates the force the second time at F2, takes a giant leap to F3, then back track again to F4. By 
reducing ti to —1/24, the back tracking steps F2 and F3 are reduced to f2 and f3. However, it is not possible to move 
any force evaluation points closer to the midpoint of the trajectory. The force evaluation points of algorithm C are 
indicated by circles. Its first and last force evaluation points nearly coincide with Fl and F4, however, it evaluates 
the force and the force gradient right at the midpoint of the trajectory as shown by the unobstructed circle. 
Better algorithms are obtained by interchanging T <^ V in 



Tm2 — ■ ■ ■ exp(eto^) exp(ewiT) exp(etiF) exp(ew27') exp(et2l^) (3.9) 
so that the momentum is updated first with the choice 



« ^ i 1 - .P^^P^ I . (3.10) 



4 V V 3ti 

Now the force is evaluated initially, at two intermediate points, and at the midpoint. The results, as shown in Fig[8] 
and Figini are much improved over the previous case. However, the locations of the forces evaluation points remain 
unusual. As shown in FiglTOl for ti — —0.1, the algorithm first evaluates the force at the starting point 2, backtracks 
past 2 to evaluate the force at F2, leaps forward to evaluate the force near the midpoint at F3, and shoots past the 
final point 3 to evaluate the force at F4. Tuning the parameter ti more negative to —0.5 reduces the back tracking 
points from F2 to f2 and F4 to f4 and improves the algorithm. However, as ti becomes even more negative, such as 
ti = — 1 or —2, those back tracking points bunch up very close to the initial and final points and do not sample the 
trajectory evenly as algorithm C. 

IV. INTEGRATORS WITH FIVE AND SIX FORCE EVALUATIONS 

For = 5, (|3.ip can be solved to give 



%{e)^... r2(aoe)T2(ai£)r2(a2£), (4.1) 

with free parameter a and coefficients 

a2 = aai, = -2^/^ {l + a^f'^ ai, (4.2) 



ai = ■ ^. (4.3) 

2(1 + a) -21/3(1 + 0,3)1/3 

The FR integrator is reproduced with a = 0. By introducing a non- vanishing 02, one is able to reduce the negative 
step size oq. This is shown in Fig[TlJ Fig [T^ and FigfT51show the energy and the precession error as a function of a. 
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While the energy error height is lowest for a = 0.5, the procession error is the smallest at a = 1. The latter is related 
to the fact that the backward step size ag is minimized at a = 1. The resulting algorithm with 

1 4I/3 
"i = ''2-j3^, ao = -j3^, (4.4) 

has long been advocated by Creutz and Gockscb^., Suzuki^ and McLachlan^°. 

Recently, a fundamental theoren*^ has allowed fourth order forward algorithms to be derived for any number of 
operators^^. In particular, one can generalize algorithm A to A'^ — 1 force plus one force-gradient evaluations (or N 
force evaluations using extrapolation). This is the class of algorithm with uniform splitting cocfhcients 

U = — - — , Vi = —7- — (4.5) 

iV-1' N{N-2y ^ ' 

and where the algorithm begins and ends with a momentum updating step: 

P' - P + (F(q) + ^^(^^^V(|F(q)p) . (4.6) 

We will denote this class of algorithm as AN. The energy and precession errors for A5 are as indicated in FigfT2l and 
FigHni Algorithm A5's precession error is more than 4 times smaller than that of algorithm C and 200 times smaller 
than non- forward algorithm (|4.ip at a = 1. As shown in FigfT4l algorithms (|4.ip again characteristically evaluates the 
force off the trajectory. As the negative time step oq is reduced by increasing a — 0.3 to a = 1.0, the off-trajectory 
force-evaluation triangle is reduced from F2-F3-F4 to f2-f3-f4. 

In FigfTSl and FigfTOl we compare the energy and precession error of A6 with that of Blanes and Moan^'^ (BM), a 
widely cited fourth-order integrator with six force evaluations. BM's energy error is comparable to that of algorithm 
C, but its maximum error height is four times that of A6. For the precession error, BM's error is two orders of 
magnitude larger than A6 and twenty times larger than C. Included in the comparison is algorithm A6', in which the 
force gradient is computed via extrapolation. Its energy error is nearly indistinguishable from that of A6, however, 
its precession error is much larger. Blanes and Moan's integrator is superior among non- forward algorithms because 
it has only two very small backward time steps, as shown in Fig ll7l (The momentum updating step first version of 
the BM integrator is not considered because it has much large errors than the position- first version discussed above.) 



V. CONCLUSIONS 



All approximation methods for solving any evolution equation should emulate its exact solution as much as possible. 
The efficiency of an algorithm cannot be decided on the basis of factorization schemes, in which only the error 
coefficient of the error commutators are known. In the past, symplectic integrators have been prized for their excellent 
conservation properties. However, because of the perceived difficulty in circumventing the Shang-Suzuki theorem, 
forward integrators were not developed until this decade. In this work, we showed that forward integrators are 
more attuned to the exact solution by evaluating the force closely on the trajectory. By comparison, non-forward 
integrators, because of their backward time steps, are constrained to evaluate the force off the trajectory, resulting in 
the loss of efficiency. In all cases studied, non-forward integrators are improved by simply reducing the size of their 
negative time steps. 
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FIG. 1: The fourth-order energy error coefficients of forward symplectic integrator A, extrapolated gradient algorithm A', non- 
forward symplectic integrator FR (Forest-Ruth) and non-symplectic integrator RKN (Runge-Kutta-Nystrom), as a function of 
time in terms of the orbital period P when solving the 2D Kepler orbit. The time step size is denoted by e. 




FIG. 2: The fourth-order orbital precession error coefRcient as measured by the rotation angle of the Laplace- Runge-Lenz 
(LRL) vector for integrators described in FiglT] 




FIG. 3: The intermediate positions and force evaluation points of the Forest-Ruth integrator. The dash curve is the exact 
orbit. See text for details. 




FIG. 4: The force evaluation points of integrator FR (solid squares), forward integrator A (circles) and non-symplectic integrator 
RKN (plus signs). 
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FIG. 5: The fourth-order energy error coefficients of a family of non-forward integrators (|3.3() with four force-evaluations 
including that of McLachlan (M) as compared to that of forward integrator C. The parameter ti characterize the size of the 
backward time step for updating the intermediate positions. 
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FIG. 6: The fourth-order orbital precession error coefficient as measured by the rotation angle of the Laplace-Runge-Lenz 
vector for integrators described in FigO 
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FIG. 7: The force evaluation points of integrator M (solid squares), reduced backward time step integrator with ti = —1/24 
(hollow squares) and forward algorithm C (circles). The solid circles denote the starting position 2 and the final position 3 
after one time step. 
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FIG. 8: The fourth-order energy error coefRcients of a family of non-forward integrators (|3.9() with four force evaluations which 
updates the momentum first. Here, the more negative the parameter ti the smaller the negative time step size for updating 
the intermediate positions. 
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FIG. 9: The fourth-order orbital precession error coefRcient as measured by the rotation angle of the Laplace- Runge-Lenz 
vector for integrators described in Fig[8] In this case, it is possible to fine tune ti so that the precession error is zero after each 
period. 




FIG. 10: The force evaluation points of integrator (|3.9|l with ti = —0.1 (sohd squares), the reduced backward time step 
integrator at ti = —0.5 (hollow squares) and the forward algorithm C (circles). 




FIG. 11: The coefficients of integrator p.lfl as a function of the free parameter a. The negative coefficient ao is least negative 
at a = 1. 




FIG. 12: The fourth-order energy error coefficients of a family of non-forward integrators (|4.1|l with five force evaluations. 
These are compared to the five-force forward integrator A5 and the four-force forward integrator C. 
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FIG. 14: The force evaluation points of non-forward integrator (|4.1|) at a = 0.3 (solid squares), with minimum backward time 
step at a = 1 (hollow squares) and that of forward algorithm A5 (solid circles). 
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FIG. 15: The fourth-order energy error coefficients of three integrators with six force-evaluations. BM is Blanes and Moan's 
integrators^. A6 and A6' are forward integrators. A6' uses the extrapolated force gradient. Algorithm C is a four-force forward 
integrator kept for comparison. 
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FIG. 16: The fourth-order orbital precession error coefficient as measured by the rotation angle of the Laplace-Runge-Lenz 
vector for integrators described in Fig llSI The extrapolated gradient integrator A6''s precession error is much larger than that 
of A6, but still smaller than that of BM. 




FIG. 17: The force evaluation points of non-forward integrator BM (solid squares) and that of forward integrator A6 (circles). 



